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ABSTRACT 

We investigate the evolution, following gas dispersal, of a star cluster produced from a hy- 
drodynamical calculation of the collapse and fragmentation of a turbulent molecular cloud. 
We find that when the gas, initially comprising m 60% of the mass, is removed, the system 
settles into a bound cluster containing ps 30 — 40% of the stellar mass surrounding by an 
expanding halo of ejected stars. The bound cluster expands from an initial radius of < 0.05 pc 
to 1 — 2 pc over w 4 — 10 Myr, depending on how quickly the gas is removed, implying that 
stellar clusters may begin with far higher stellar densities than usually assumed. With rapid 
gas dispersal the most massive stars are found to be mass segregated for the first ~ 1 Myr 
of evolution, but classical mass segregation only develops for cases with long gas removal 
timescales. Eventually, many of the most massive stars are expelled from the bound cluster. 
Despite the high initial stellar density and the extensive dynamical evolution of the system, 
we find that the stellar multiplicity is almost constant during the 10 Myr of evolution. This is 
because the primordial multiple systems are formed in a clustered environment and, thus, by 
their nature are already resistant to further evolution. The majority of multiple system evolu- 
tion is confined to the decay of high-order systems (particularly quadruple systems) and the 
formation of a significant population of very wide (10 4 — 10 5 AU) multiple systems in the 
expanding halo. This formation mechanism for wide binaries potentially solves the problem 
of how most stars apparently form in clusters and yet a substantial population of wide binaries 
exist in the field. We also find that many of these wide binaries and the binaries produced by 
the decay of high-order multiple systems have unequal mass components, potentially solving 
the problem that hydrodynamical simulations of star formation are found to under-produce 
unequal-mass solar-type binaries. 

Key words: binaries: general - methods: N-body simulations - stars: formation - stars: low- 
mass, brown dwarfs. 



1 INTRODUCTION 

For many decades, it has been possible to perform A'-body simula- 
tions of the evolution of stellar clusters, examining their dynamical 
evolution and the evolution of their stellar populations ( [von Ho-] 
|erner|1960|[van Albada|1968l|Aarseth & Woolf|1972^ . Some early 
work was even done on the dynamical evolution of young clus- 
ters. For example, Aarseth & Hills (1972) studied the dynamical 
evolution of a stellar cluster with initial sub-clustering and found 
that sub-clustering increased both the fraction of stars ejected dur- 
ing the relaxation of the cluster and the formation of binary stars. 
However, in the last two decades, during which observational stud- 
ies of star-forming regions have dramatically increased in number, 
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much more attention has been paid to the dynamical evolution of 
young stellar clusters. 

|Kroupa| ( |1995a} hypothesised that all stars were formed in 
binary systems and studied whether such initial conditions could 
lead to the observed field population of binaries. With an initial 
period, P, distribut ion that was flat in log P ranging from 10 3 to 
10 7 ' 5 days, |Kroupa| found that subsequent dynamical interactions in 
dense stellar clusters could reduce the binary fraction from 100% 
and modified the period distribution to match the properties found 
for field stars. On the other hand, simulations that begin only with 
binaries cannot reproduce the observed frequencies or properties of 
triples ( Portegies Zwart et al. 2004 1 and simulations that begin with 
primordial binaries and triples cannot reproduce the observed fre- 
quencies of quadruples ( |van den Berk, Port egies Zwart & McMil- 
|lan|2007) . From his simulations of modest-sized clusters contain- 
ing 200 binaries, Kroupa ( 1995b) also found that the clusters would 
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be surrounded by a halo of ejected stars with a lower binary frac- 
tion. |Kroupa (1998) found that ejected binary systems tended to 
be massive with a preference for equal-mass components. He also 
found that the mean system mass was approximately independent 
of ejection velocity from 2-30 km s~\ with dynamical ejection 
providing relatively massive stars as far as 300 pc from the cluster 
within 10 Myr (see also Sterzik & Durisen 1995 ). The binary frac- 
tion decreased monotonically with ejection velocity for low-density 
clusters, but for dense clusters the dependence was more complex. 
However, long-period systems (periods > 10 4 days) could not be 
ejected with large velocities (> 30 km s _1 ). 

More recently, attention has turned to the effects of gas disper- 
sal on the evolution of young clusters and, following the pioneer- 
ing work of Aarseth & Hills ( 1972), to the effects of sub-structure. 
Simple analytical models indicate that a young cluster requires a 
star formation efficiency of 77 > 50% to survive the rapid dis- 
persal of gas as a bound cluster (Hills 1980). However, more so- 
phisticated semi-analytic and numerical models show that the sur- 
vival of a bound cluster is more complex. If the gas is removed 
on a timescale significantly longer than the crossing time, then r\ 



as low as 20% can still result in a bound cluster (Elmegreen 1983 



Mathi eu|1983||Verschueren & David|19 89 ; Geyer & Burker t|2001| 
Baumgardt & Kroupa 2007 1 and even if the gas is removed quickly 



many stars can be lost, but a bound core can remai n (|Lad a, Margulis 
& Dearbom|1984||Goodwin|1997||Adams|2000||Geyer & Burkert| 
2001||Kroupa, Aarseth & Hurley|2001||Boiry & Kroupa|2003a|b| 
Baumgardt & Kroupa 2007 1. For example, Geyer & Burkert ( 2001 1 

0.4 and rapid ex- 



Boily & Kroupa (2003b ) find that with 77 
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plusion around 30-40% of the initial cluster could remain bound. 
|Baumgardt & Kr oupa (2007) find that for instantaneous gas re- 
moval r\ > 1/3 is required for such a remnant cluster to be formed, 
but that if the gas is removed over several crossing times, around 
50% of the cluster can remain bound even for r\ as low as 15-20%. 
Finally, Goodwin ( 2009 1 notes that the survivability of a cluster de- 
pends not only on the star formation efficiency, but on the cluster's 
dynamical state. In particular, if the gas removal occurs when the 
cluster is still dynamically cool (subvirial), then cluster survivabil- 
ity can be greatly increased. 

The dynamical state of a cluster is also important for mass 
segregation. Bonnell & Davies (1998 ) showed that the Orion Neb- 
ula Cluster (ONC) is not old enough for the observed segregation 
of its massive stars to be due to dynamical segregation alone. In- 
stead, they argued that primordial segregation was required. Re- 
cently, however, McMillan ~et al.H2007) , |Fellhauer et al.| ( |2009| >, 
and [Moeckel & Bonnell (2009) have shown that if the cluster is 
formed via the merging of sub-clusters, mass segregation can be 
accelerated. Similarly, Allison et al. (2009 ) investigated the evolu- 
tion of dynamically cool clusters with fractal initial structure and 
found that dynamical mass segregation of the most massive stars 
can occur quickly. 

However, the above studies have begun with artificial initial 
conditions. Assumptions must be made about the initial structure of 
the stellar clusters (e.g. Plummer spheres, fractal or sub-clustered 
stellar distributions), the virial ratio (e.g. in virial equilibrium, sub- 
virial so they collapse, or super-virial so they expand), the initial 
populations of binary (e.g. 100% binaries) and higher-order multi- 
ple systems (typically, no triples or higher-order systems, with the 
exception of |van den Berk et al.|2007 ), and the properties of these 
multiples such as their distributions of separations and mass ratios. 

Recently, it has become possible to perform hydrodynamical 
simulations of the collapse of a molecular cloud to form a star clus- 
ter containing in excess of 1000 of stars and brown dwarfs while re- 



solving the opacity limit for fragmentation (Bate 2009a). Such cal- 
culations allow detailed comparisons to be made with the observed 
properties of young clusters, such as the stellar initial mass function 
(IMF), multiplicity, the properties of binary and multiple systems, 
and the global properties of stellar clusters. Bate's calculation in- 
cluded gravity and hydrodynamics using a barotropic equation of 
state and used sink particles to model the stars and brown dwarfs. 
The sink particles had accretion radii of only 5 AU, allowing cir- 
cumstellar discs to be resolved down to ~ 10AU in radius, while 
binaries and multiple systems with separations as close as 1 AU 
could be followed. Despite the limited amount of physics included 
in the calculations, many of the stellar properties obtained from the 
calculation were in good agreement with the observed properties 
of stellar systems. For example, the observed stellar multiplicity as 
a function of primary mass was found to be well reproduced, as 
were the separations and mass ratios of low-mass stars and brown 
dwarfs. Even the distribution of the orientation of the orbital planes 
of triple stellar systems was found to be in good agreement with ob- 
servations. The two main deficiencies of the calculation were that 
it produced a ratio of brown dwarfs to stars well in excess of that 
which is observed, and solar-type binaries were found to have a 
deficit of low mass ratio systems. The first of these is likely to be 
due to the neglect of radiative feedback from the protostars piate| 
|2009b||Offher et al.|2009) , while the solution to the latter problem 
is currently unclear. However, the implication of Bate's calcula- 
tion is that many stellar properties result primarily from dissipative 
gravitational dynamics and that additional physical processes such 
as radiative transfer and magnetic fields may not be very important 
for determining these properties. 

One objection that could be raised to Bate ( 2009a) is that, due 
to computational limitations, the formation of the cluster could only 
be followed for 2.85 x 10 5 yr (1.5 global free-fall times) which 
was only 1.5 x 10 yr after the first star formed. However, in many 
cases, the stellar properties were compared to star-forming regions 
that were typically ~ 1 Myr old, or even to field stars. When Bate's 
calculation was stopped, there were many triple, quadruple, and 
higher-order stellar systems, and many of there systems were dy- 
namically unstable. Furthermore, of the 500 Mq of gas originally 
contained in the molecular cloud only 191 Mq (38%) had been 
converted to stars by the end of the calculation. Star formation 
would have continued if the calculation were continued. But even 
if the gas was dispersed at this point and star formation ceased, 
there is the question of how the stellar properties would evolve over 
million-year timescales and longer. 

In this paper, we take the end point of the hydrodynamical 
simulation of Bate (2009a) as the starting point for yV-body simula- 
tions. The advantage of these starting conditions is that the stellar 
cluster has been formed self-consistently from the evolution of a 
gravitationally-unstable turbulent molecular cloud. Thus, the struc- 
ture of the cluster and the properties of binary and multiple systems 
have all been consistently set up, albeit with a restricted number of 
physical processes (e.g. no radiative feedback or magnetic fields). 
Our aims are to determine the long-term (10 Myr) evolution of the 
cluster and how its properties and those of the stellar multiple sys- 
tems evolve during the dispersal of the gas in the cluster. The only 
previous similar work we are aware of is that of Hu rley & Bekki| 
(2008 ) who also performed SPH simulations of star formation to 
generate the initial conditions for yV-body simulations. However, 
their work studied star cluster formation on giant molecular cloud 
scales, assumed equal-mass stars, and did not resolve primordial 
multiple systems. 

The paper is structured as follows. Section 2 briefly describes 
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the initial conditions and the numerical method for the calculations. 
The results are discussed in Section 3. In Section 4, we discuss 
the implications of the results for our understanding of star cluster 
evolution. Our conclusions are given in Section 5. 



2 COMPUTATIONAL METHOD 

The calculations discussed in this paper are continuations of a hy- 
drodynamical calculation that began with a turbulent molecular 
cloud and produced a star cluster. The hydrodynamical calculation 
was performed using a smoothed particle hydrodynamics (SPH) 
code (Bate, Bonnell & Price 1995) and the cloud was evolved for 
1.5 initial cloud free-fall times (2.85 x 10 yr, where the free-fall 
time was tg = 6.0 x 10 12 s or 1.90 x 10 5 yr). The stars and brown 
dwarfs were modelled using sink particles. In this paper, we take 
the sink particles at the end point of the hydrodynamical calculation 
as the starting point for a series of A'-body simulations and evolve 
the stellar cluster to an age of 10 Myr. Note that for all the times we 
quote, the starting point is taken to be that of the hydrodynamical 
calculation, not the A'-body simulations, unless stated otherwise. 



2.1 Initial conditions 

The initial conditions for the A'-body simulations are taken from 
the main star cluster formation calculation of Bate (2009a|. The 
initial conditions for the hydrodynamical calculation consisted of 
a 500 Mq, uniform, spherical molecular cloud represented by 
3.5 x 10 7 SPH particles. The cloud's radius was 0.404 pc (83 
300 au). At the initial temperature of 10 K, the mean thermal 
Jeans mass was 1 Mq (i.e. the cloud contained 500 thermal Jeans 
masses). Although the cloud was uniform in density, an initial su- 
personic turbulent velocity field was imposed in the same man- 
ner as |Ostriker, Stone & Ga mmie (2001 1. A divergence-free ran- 
dom Gaussian velocity field was generated with a power spectrum 
P(k) oc fe _4 ,where k is the wavenumber. In three dimensions, 
this results in a velocity dispersion that varies with distance, A,as 
(j(A) oc A 1 / 2 , in agreement with the observed Larson scaling re- 
lations for molecular clouds jLarson|198 1 ). The velocity field was 
normalized so that the kinetic energy of the turbulence equaled the 
magnitude of the gravitational potential energy of the cloud. Thus, 
the initial rms Mach number of the turbulence was M — 13.7. The 
'turbulence' decayed as the calculation was evolved; there was no 
'turbulent driving'. The cloud was allowed to evolve freely into the 
vacuum surrounding it; there were no boundary conditions applied 
to the simulation. 



2.1.1 Sink particles 

To model the thermal behaviour of the gas without performing ra- 
diative transfer, the SPH calculation used a barotropic equation 
of state (see |Bate et aL] {2003) for further details). This equa- 
tion of state kept the temperature at 10 K for densities < 10~ 13 
g cm~ 3 , while above this density the temperature increased with 
density in order to mimic the opacity limit for frag mentation |Low| 
& Lynden-Beli] [19761 |Rees][lT>76l [Sllk||1977a|b[ |Boyd & Whit? 



worth||2005^ . When the central density of a fragment exc eeded 
„ in- 11 ~ — , — 3 - t i„ rt — i u.. „ „J«1, 177777 
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gem , it was replaced by a sink particle iBate 



|etal.|1995| l. In the main calculation of Bate ( 2009a), sink particles 
were formed by replacing the SPH gas particles contained within 
r^cc = 5 AU of the densest gas particle in a pressure-supported 
fragment by a point mass with the same mass and momentum. Any 



gas that later fell within this radius was accreted by the point mass 
if it was bound and its specific angular momentum was less than 
that required to form a circular orbit at radius r acc from the sink 
particle. Sink particles interacted with the gas only via gravity and 
accretion. The gravitational acceleration between two sink parti- 
cles was Newtonian for r ^ 4 AU, but was softened within this 
radius using spline softening Benz ( 1990). The maximum acceler- 
ation occurs at a distance of « 1 AU; therefore, this is the minimum 
separation that a binary had in the hydrodynamical simulation even 
if, in reality, a binary's orbit would have been hardened. Sink par- 
ticles were permitted to merge if they passed within 2Rq of each 
other; only one merger occurred during the simulation. 

2.1.2 The end state of the hydrodynamical calculation 

The end state of the hydrodynamical calculation consisted of a 
cluster of 1253 stars and brown dwarfs with a combined mass of 
191 M while 309 Mq of gas remained. Note that [Batel p009"al > 
states that there were 1254 objects, but in fact it was not noticed 
when the paper was written that two of these had merged and had 
been replaced by a single object. This has no significant effect on 
the statistics that were reported in the original paper. The final clus- 
ter was created from the merger of around half a dozen sub-clusters 
that formed in gaseous filaments that were produced in the turbu- 
lent gas. The sub-clusters formed in a dynamically cool configu- 
ration and fell together to produce a single cluster surrounded by 
an extended halo of unbound, dynamically ejected stars and brown 
dwarfs. This cluster was very dense, with a half-mass radius of only 
10 4 AU (0.05 pc) and a velocity dispersion of around 4 km s _1 , 
with the velocities of the ejected halo objects extending beyond 
20 km s -1 . Many of the stars were in binary or higher-order multi- 
ple systems and the stellar multiplicity was found to be a strongly 
increasing function of primary mass, with values in good agree- 
ment with observed values. Many of the trends for the separation 
and mass ratio distributions of observed binaries were also repro- 
duced (e.g. the trend for very-low-mass binaries to have smaller 
separations and mass ratios nearer unity than for low-mass and 
solar-type stars), with the exception that unequal-mass solar-type 
binaries were under-produced. 

Further details regarding the initial conditions, the SPH code, 
the evolution of the hydrodynamical star cluster formation calcula- 
tion and the statistical properties of the stars and brown dwarfs can 
be found in Bate ( 2009a). We now turn to the question of how this 
initial cluster evolves on 10 Myr timescales. 

2.2 The A'-body calculations 

The final positions, masses and velocities of the sink particles in 
the hydrodynamic simulation are taken as the initial conditions for 
a gravitational A'-body simulation. The integration was performed 
with the NBODY6 code (Aarseth 2003), which uses a fourth-order 
Hermite integrator. We advanced the system for 10 Myr from the 
end of the hydrodynamic simulation, and the relative energy error 
was of order 10~ 5 . It is usual practice when performing A'-body 
calculations with a modest number of stars to generate many sets of 
initial conditions and consider ensemble quantities, so that general 
trends can be separated from chaotic variations in the integrations. 
In this case, with initial conditions generated not by random real- 
izations of assumed bulk cluster properties but by a self-consistent 
and expensive hydrodynamical calculation, we are forced to take 
a slightly different approach to disentangling general results from 
those particular to the exact initial conditions. 
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Calculation 


End Time 


Remnant Cluster 


Mass 




Final Numbers of Syst 


ems 




[Myr] 


Mass 


Radius [pc] 


Segregation 


Single 


Binary 


Triple 


Quadruple 


Hydrodynamical (Bate 2009a) 


0.3 






10 most masssive stars 


905 


90 


23 


25 


Instant Gas Removal 


10 


30% 


2 




956 


98 


23 


8 


T .i = 0.3 Myr 


10 


30% 


1.8 


J> 1 Mg 


941 


100 


20 


13 


T .i = lMyr 


10 


40% 


1.3 


^1M 


962 


97 


19 


10 


No Gas Removal 


10 


50% 


0.8 


^ 0.3 M 


966 


93 


11 


17 



Table 1. A summary of the calculations discussed in this paper. The four W-body calculations with different gas removal timescales take their initial conditions 
as the stars (sink particles) at end of the hydrodynamical calculation of Bate 1 2009a). Tq.i is the time from the start of the JV-body evolution at which the gas 
mass has fallen to 10% of its initial value. In addition to the calculations listed above, the instantaneous gas removal case was run ten times with randomly 
perturbed initial conditions. In the table, we summarise the properties of the remnant clusters that are left at the end (10 Myr) of the W-body calculations (their 
radii and the mass fractions of the stars they contain). We also comment on the mass segregation that is observed (e.g. stars with masses > 0.3 Mq at the end 
of the case with no gas removal). Finally, we give the numbers of single and multiple systems at the end of the hydrodynamical calculation of Bate 1 2009a I 
and each of the unperturbed W-body calculations. The W-body evolution always results in a significant decrease in the number of quadruple systems and small 
increases in the numbers of single and binary systems. Slower gas removal timescales generally produce more compact and populous remnant clusters that 
show more mass segregation. 
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Figure 1. The solid line gives the cumulative radial distribution of the 
309 Mq of gas that remained at the end of the hydrodynamical calculation 
of Bate (2009a), although the gas distribution is not spherically symmetric. 
Note that the half-mass radius of the stellar cluster is only 0.05 pc. The dot- 
ted line gives the cumulative radial distribution of a Plummer sphere with 
the same total mass and a characteristic radius of 0.25 pc. The two distribu- 
tions have similar half-mass radii. 

Ten sets of perturbed initial conditions were generated by tak- 
ing the initial conditions and randomising each body's position by 
10 -4 times its radius from the cluster's center of mass. That is, each 
position vector r was adjust so that r — > r + 10~ 4 rw, with w a 
randomly oriented unit vector. The initial specific potential energy 
of each body was thus perturbed by an amount greater than the inte- 
gration accuracy. Stars in binaries were treated as single bodies for 
the perturbation, so that their initial binary properties remained the 
same. By comparing results of the main integration to those of the 
perturbed versions of the initial conditions, we can ascertain which 
long-term results are robust and general to the type of cluster gen- 
erated by the hydrodynamic simulation, and which are sensitive to 
the exact state of the cluster. 

The main calculation and the ten randomly perturbed calcu- 
lations are all run assuming that the gas left at the end of the hy- 
drodynamical calculation is instantaneously removed. This is the 
most extreme gas removal history possible and, based on previ- 
ous work, should result in the least-bound cluster remnant possible 
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Kroupa 2007). If the gas is removed over a finite period of time, 
we would expect more of the cluster to remain bound and it might 
also affect the frequencies and properties of the multiple stellar sys- 
tems. To investigate how much of a difference slow removal of the 
gas makes, we have performed three additional calculations that in- 
clude an additional gravitational potential to model the contribution 
of the gas. 

The gas mass at the end of the hydrodynamic calculation is 
309 Mq. While the actual gas distribution is not spherical (see Fig- 
ure 1 of Bate 2009a), we approximate it as a Plummer sphere with 
scale radius 0.25 pc so that the half-mass radius is similar to that 
of the actual gas distribution (see Figure[T|(. The mass of the Plum- 
mer potential is decreased with time, t, from the beginning of the 
A-body calculations so that the gas mass is given by 



where A/o is the original gas mass (309 Mq) and T ga3 is a gas 
removal timescale. We ran three cases: moderately fast gas disper- 
sal, where the mass was reduced to 10% of its initial value after 
0.3 Myr of A-body evolution (total age 0.6 Myr); slow dispersal, 
where the mass was reduced to 10% at 1 Myr of A-body evolution 
(total age 1.3 Myr); and an (unphysical) test case where the gas 
mass remained constant throughout the calculation. 

A summary of the main properties of the calculations dis- 
cussed in the next section is presented in Table [T] The plots and 
discussion that follow are based on the unperturbed initial condi- 
tions with instantaneous gas removal, with comparison to the ran- 
domized integrations and those with gas potentials made where ap- 
propriate. 



3 RESULTS 

3.1 Evolution of cluster properties 

3.1.1 Substructure 

In the hydrodynamic simulation of Bate (2009al, the stars form in 
a structured fashion, with smaller subclusters merging to form a fi- 
nal cluster consisting of a tightly bound core with radius ~ 0.05 pc 
surrounded by an expanding halo of ejected stars. The reader inter- 
ested in the hydrodynamical evolution of the cluster is encouraged 
to consult Figure 1 of [Bate]p009a) . In Figure [2] of this paper, we 
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Figure 2. The iV-body evolution of the unperturbed star cluster with instantaneous gas removal. Each star is plotted using an open circle with the area of the 
circle proportional to the star's mass. Stars with masses ^ 3 Mq are plotted in red, while those with masses 1 ^ M < 3 Mq are plotted in blue. Binaries 
and multiples can be seen as almost concentric circles when the components have different masses. We show the structure on 200, 20, and 5 pc scales at four 
different times. On the largest scales, it can be seen that some ejected stars in the halo reach distances > 100 pc in less than 10 Myr and that these stars are 
a mixture of massive and low-mass stars. On the smallest scales, the remnant bound cluster containing approximately 1/3 of the stars expands from its initial 
radius of si 0.05 pc to pb 2 pc. The most massive star and a small cluster of 20 companion stars is ejected from the remnant cluster at ss 1 Myr and can be 
seen to the lower-left of the main cluster remnant in the lower two 5 pc panels. 
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Figure 3. The value of the parameter Q over the first 1 .5 Myr. At each time 
the quantity is determined for three orthogonal projections of the cluster, 
shown as the points, and the means of those values are connected by the 
lines. The black line gives the evolution of the entire stellar cluster and 
halo, while the grey line gives the evolution of the inner 3 pc region. The 
inner region evolves from a sub-clustered stellar distribution (Q < 0.8) to 
a core and halo structure (Q > 1) during the hydrodynamical calculation. 
When the gas is removed, the remnant core expands to fill most of the inner 
3-pc and Q decreases to approximately unity. 

display the global evolution of the stellar cluster during the 10 Myr 
of A'-body evolution. Each star is represented by an open circle with 
its area proportional to the star's mass. On large scales, the stars in 
the halo are seen to expand to distances exceeding 100 pc from 
the cluster in less than 10 Myr. On small scales, a bound cluster 
remnant containing approximately 30% of the stars expands from 
a radius of ~ 0.05 pc to ~ 2 pc over several Myr. 

The degree of subclustering can be quantified by the parameter 
Q qCartwright & Whitworth||2004| >. Values of Q < 0.8 indicate 
substructure and can be associated with a fractal dimension, while 
values greater than 0.8 are associated with radial density variation. 
Calculating Q requires both an area and radius to be assigned to a 
cluster; as the cluster at early times is decidedly non-spherical, we 
calculate the area as the convex hull of the stars' positions projected 
onto a plane, and take the radius to be that of a circle having the 
same area. The evolution of Q over the first 1.5 Myr of the cluster 
evolution is shown in Figure [5] 

At the earliest measured times the cluster consists of 154 stars, 
and the value of Q ~ 0.4 corresponds to an estimated fractal di- 
mension of 1.5, i.e. highly substructured. Over the next ~4x 10 4 
yr the early subclusters merge, the number of stars increases to 579, 
and Q rises above 0.8. This indicates a fractal dimension of 3, i.e. a 
distribution lacking any significant substructure. By the end of the 
hydrodynamical calculation the number of stars has reached its fi- 
nal value of 1253, and Q has risen in excess of 1.5. In |Cartwright| 
& Whitworth (2004 ), larger values of Q are identified as corre- 
sponding to an increasingly steep radial density profile of a spheri- 
cal cluster. In our case the cluster is not well-fit by a single power 
law density profile, but rather has a core-halo structure, which com- 
plicates a straightforward interpretation of Q. Because the cluster 
is characterized by a bound remnant surrounded by an expanding 
halo, the exact value we calculate is quite sensitive to the radial ex- 
tent of stars that we consider. In Figure [3] following gas removal, 
two lines are plotted. The black line includes all the stars, arranged 
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Figure 4. Lagrangian radii during the span of the main A'-body calculation. 
From bottom to top, solid curves enclose 0.025, 0.05, 0.1. 0.2, 0.3, 0.4, 
0.5, 0.7, and 0.9 of the total stellar mass. The dotted line shows the median 
radius of the five most massive stars. 

with a core-halo structure, while for the grey line the analysis has 
been restricted to the central 3-pc of the cluster. In this central re- 
gion, the cluster core expands upon gas removal to occupy most 
of the inner region and the value of Q drops to unity, reflecting an 
absence of significant structure and the fact that the halo is only 
present on larger scales. For our purposes it suffices to note that Q 
remains quite high for the entire 10 Myr run, showing the persis- 
tence of the core-halo structure. 

We emphasise that the evolution of the cluster from a highly- 
clustered to a core-halo structure occurred while the the cluster was 
still embedded and forming; the number of stars increases by nearly 
an order of magnitude during the transition from a fractal dimen- 
sion ~ 1.5 to ~ 3.0. This highlights the importance of coupled 
gas and gravitational dynamics in young cluster evolution, and the 
difficulty in assigning realistic initial conditions to purely A'-body 
studies. 

3.1.2 Lagrangian radii and stellar velocities 

The radial mass structure of the cluster can be seen in the evolution 
of the Lagrangian radii, that is the radii encompassing a constant 
percentage of the cluster's mass. This calculation is dependent on 
what one considers the cluster center. At each time we identify the 
star which is located in the region of the greatest stellar density, 
which we define as the star with the smallest sphere encompassing 
50 neighbours, and set this star as the origin. The behavior of these 
diagnostics depends sensitively on the details of the integration, 
especially the behavior of the most massive stars. We first discuss 
the main features of the Lagrangian radii for the unperturbed initial 
conditions, shown in Figure]?] and then note how these vary among 
the randomly perturbed clusters and those with a gas potential. 

The first Myr of A'-body evolution is marked by general ex- 
pansion of the outer 70% of the cluster's mass, with the Lagrangian 
radii increasing by an order of magnitude. In contrast, the inner 
10% of the mass expands by only a factor of 2. Between 1 and 
2 Myr these inner radii expand by an order of magnitude, catch- 
ing up to the expansion of the outer radii. At later times there is 
a divergence as the outer radii continue to expand, while the inner 
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Figure 5. The evolution of the stellar velocity dispersion and the radius of 
the bound remnant cluster in the main iV-body calculation. In the top panel 
we give the time evolution of the radius containing 30% of the stars (dashed 
line), centred on the star at the location of the highest stellar density. It can 
be seen that the remnant cluster expands by a factor of fti 20 in size. The 
velocity dispersion of the stars in the remnant (taking the centre of mass 
velocity for all binaries and counting them only once) is given by the solid 
lines in both panels. The velocity dispersion decreases from ss 4 to ss 1 
km -1 . The spike at 1 Myr is associated with the ejection of the most 
massive binary from the cluster. In the lower panel, the velocity dispersion 
is plotted as a function of the remnant cluster's radius. The dotted line has 
a slope of — 1 /2, as expected if the remnant cluster remains close to virial 
equilibrium as it expands. 



radii (certainly the inner 20%, and after 8 Myr the inner 30%) halt 
their growth, leaving a remnant cluster core at the center of an ex- 
panding halo. The evolution of this remnant cluster core is further 
explored in Figure [5] where we plot the velocity dispersion of the 
stars in the remnant as a function of both time and of the remnant 
cluster's radius. As expected for a system near virial equilibrium, 
the velocities of the stars decrease as the remnant expands roughly 
as iT 1 / 2 . 

The 1 Myr delay in the expansion of the inner 10% of the mass 
is explained by the behavior of the most massive stars. The 5 most 
massive stars in the cluster are by themselves approximately 10% 
of the total stellar mass, and if they are concentrated in the clus- 
ter center then they will dominate the inner Lagrangian radii. To 
give an idea of these stars' radii from the cluster centre, we also 
plot in Figure [4] the median radius of the 5 most massive stars as 
a dotted line. At about 1 Myr a massive binary, comprised of a 
5.35 Mq primary (also the most massive star in the cluster) and a 
3.54 Mq secondary, is ejected from the cluster center (see the lower 
two right-hand panels of Figure [2} . This event also appears in the 
velocity dispersion evolution depicted in Figure[5] The massive bi- 
nary is accompanied by a small group of ~ 20 companion stars. 
The inner 10% Lagrangian radii track this binary and its compan- 
ions for a fewx 10 4 years, at which point they begin to reflect the 
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Figure 6. Lagrangian radii during the span of the W-body calculation that 
includes a gas potential which decreases to 10% of its initial value over 
1 Myr of evolution (total age 1.3 Myr). From bottom to top, solid curves 
enclose 0.025, 0.05, 0.1. 0.2, 0.3, 0.4, 0.5, 0.7, and 0.9 of the total stellar 
mass. The dotted line shows the median radius of the five most massive 
stars. In comparison with Figure|4] the remnant cluster is smaller and more 
populous, while more of the most massive stars have been ejected. 



overall cluster structure rather than the positions of a few massive 
stars. 

This feature, a delayed expansion of the inner Lagrangian radii 
due to a clustering of the most massive stars, is seen in 9 of the 10 
randomly perturbed simulations, though the time of the ejection 
ranges from ~ 4 x 10 4 yr to 2 Myr after the start of the W-body 
evolution. In one case the massive stars remain in the cluster cen- 
ter throughout the 10 Myr. In other words, in only one of our 11 
runs with instantaneous gas removal do the massive stars remain 
continually associated with the center of the cluster for more than 
2 Myr. In two cases the massive stars return to the center of the 
cluster, one of which sees a second ejection. In that perturbed re- 
alization, after the second ejection of massive stars at about 3 Myr 
the inner Lagrange radii expand freely. While a cluster remnant is 
still identifiable, it is globally dispersing unlike any of the other 
realizations. Another feature is the behavior of the 20% radius in 
the unperturbed case, which is fairly flat after its initial expansion. 
In six of the perturbed runs this radius increases by a factor of 3 
or less between 2 and 10 Myr; in two cases this radius is flat for 
most of the evolution but grows over the final 2 Myr; and in one 
case it expands freely after the ejection of the massive stars. The 
global evolution of the cluster's mass profile thus appears to be sur- 
prisingly sensitive to the initial conditions. We can say with some 
confidence that the massive stars are unlikely to remain associated 
with the main cluster remnant, but the final characteristics of that 
remnant (i.e. its radius, mass, and its stability or expansion) display 
significant variation. 

As expected, the A'-body calculations that include a gas po- 
tential and slower gas removal result in smaller final Lagrangian 
radii, and larger membership of the remnant cluster. In Figure [6] 
we show the evolution of the Lagrangian radii in the slow gas dis- 
persal case (where the gas mass is reduced to 10% of its initial 
value after 1 Myr of iV-body evolution). Comparing this with Fig- 
ure [4] we find that somewhat more of the stellar mass (about 40% 
compared to 30%) remains in a remnant cluster and this remnant 
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Figure 7. Lagrangian radii at the end of the iV-body calculations (10 Myr) 
as a function of the gas removal timescale, To. i , which gives the time from 
the beginning of the A'-body evolution at which the gas potential has been 
reduced to 10% of its initial value (0, 0.3, 1 Myr and no gas removal). 
Longer gas removal times yield more compact and populous clusters. 

cluster is somewhat smaller in the slow gas dispersal case (radius 
approximation 1.3 pc rather than 2 pc). The effect of the gas po- 
tential is more clearly demonstrated in Figure [7] where we show 
the Lagrangian radii at 10 Myr for each of the unperturbed cases, 
which we label by the timescale, To.i, for the gas mass to be re- 
duced to 10% of its initial value. The inner Lagrangian radii all 
display the expected trend, with longer gas-removal times yielding 
a more compact cluster and populous cluster. At 10 Myr the instan- 
taneous, 0.3 Myr, 1.0 Myr, and no gas removal cases have 552, 595, 
756, and 944 stars within 3 pc of the densest star, respectively. 

3.2 Mass segregation 

Mass segregation is seen in several young clusters, typically ap- 
pearing as a concentration of the most massive stars (e.g. the ONC 
or Mon R2 |Hillenbrand & Hartmann|1998[|Carpenter et al.|1997) 
rather than a general trend affecting all mass ranges, as would 
be expected in older, dynamically mature clusters. Recent work 
(McMillan et al.|2007| [Fellhauer et al.l2009"l |Ailison et al.|2009| 
|Moeckel & Bonn ell 20091 shows that a structured mode of star 
formation, with merging substructures forming the final quasi- 
spherical cluster, can lead to segregation affecting only the massive 
stars, mirroring observations. Moeckel & Bonnell (2009 ) noted that 
the final state of the Bate (2009a) hydrodynamic simulation showed 
a similar mass segregation signal to the ONC. Here, we examine 
how this rapid mass segregation evolves on longer timescales. 

With the realization that mass segregation in very young clus- 
ters is qualitatively different from the segregation that develops as 
a spherical cluster evolves over many relaxation times, new meth- 
ods of quantifying segregation have been introduced. We use two 
diagnostics here: first, a standard technique, the cumulative radial 
distribution of stars in different mass bins; second, a variation of 
the minimum spanning tree (MST) method introduced in | Allison] 
|et al.|p009] >, in which we measure the effective surface density of 
the most massive stars and compare it to the expected surface den- 
sity. This is the method used in Moeckel & Bonnell (2009}, and 
full details are found there. Briefly, the expected surface density of 
n stars in the cluster is found by calculating the area of the con- 



vex hull, H c (n), for many samples of n random stars. The median, 
±lcr, ±2<7, and ±3<r of the quantity H c (n) jn define the 'normal' 
surface density of n random stars. This distribution is then com- 
pared to the surface density of the n most massive stars. With both 
of these methods, we restrict our analysis to stars within 3 pc of 
the star at the location of greatest stellar density. This radius en- 
compasses the final cluster remnant at 10 Myr, and prevents very 
distant escaping stars from interfering with a sensible interpretation 
of the results. 

In Figure [8] we show the state of the cluster's segregation at 
three times; the beginning of the A'-body simulation, just after 1 
Myr when the most massive stars are about to be ejected from the 
cluster center, and at around 2 Myr, which is representative of later 
times. At the earliest time, the cumulative distribution shows very 
little evidence of mass segregation (as noted by Bate 2009a), while 
the surface density analysis shows that the 10 most massive stars 
have a surface density higher than the 3a value expected for the 
cluster. At around 1 Myr, the cumulative distribution shows a clear 
difference in the radial distribution of stars above 1 Mq. The sur- 
face density of the four most massive stars is still very high, and 
both plots reflect the lowered surface density of the cluster after its 
initial expansion. It is interesting to note that the surface density 
of the cluster between 1 and 2 Myr has dropped by only a factor 
of ~ 2, while the inner Lagrangian radii (Figure]?} have increased 
by over an order of magnitude, emphasizing that the most massive 
stars were dominating the inner Lagrangian radii at early times. 

By 3 Myr, mass segregation among the four most massive stars 
remaining in the cluster has reestablished itself in the surface den- 
sity diagnostic, and remains for the remainder of the simulation. 
However, in contrast to the segregation at 1 Myr, this is due to the 
presence of two massive and well-separated binaries, and not a real 
clustering of four stars. The disappearance and return of the mass 
segregation signal is dependent on what one considers part of the 
cluster versus the field. In this case, when dispersing massive stars 
cease to be considered part of the cluster (by our 3 pc selection cri- 
terion), the remaining massive stars are left in a configuration that 
appear to be clustered. In addition to the cluster membership deter- 
mination, the amount of segregation seen in a cluster depends on 
the viewing angle. A projection that minimizes the separation be- 
tween two massive binaries will overemphasize the clustering. This 
is important to bear in mind when considering clustering involving 
only a handful of massive stars. 

In Table [2] we examine the radial properties of the remnant 
cluster and halo of ejected stars using logarithmically-spaced bins. 
In agreement with the diagnostics discussed above, there is no ev- 
idence of mass segregation. Some of the most massive stars in the 
system are found in each radial bin. The velocities of the ejected 
stars in the halo increase with radius outside of 3 pc from the clus- 
ter centre (i.e. faster moving stars have travelled further from the 
cluster centre). 

All four mass-removal scenarios begin with the same initial 
stellar positions and, therefore, all begin with the 10 most mas- 
sive stars being centrally concentrated. While this primordial seg- 
regation is transient and only to be observed in the earliest stages 
of cluster evolution, the more compact remnant clusters that result 
from slower gas removal have a direct effect on the later mass seg- 
regation. Calculated at 10 Myr, the half-mass relaxation times of 
the calculations with instantaneous gas removal, To.i = 0.3 Myr, 
To.i = 1 Myr, and no gas removal are 20.8 Myr, 12.3 Myr, 8.2 Myr, 
and 3.3 Myr, respectively. Thus, whereas the main calculation is 
dynamically unrelaxed, the two calculations with slow or no gas 
removal have evolved for more than one relaxation time, which is 
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Figure 8. The mass segregation in the unperturbed cluster with instantaneous gas removal at three times; the end state of the hydrodynamic simulation, just 
before the ejection of the most massive stars from the cluster center, and at ~ 2 Myr. Top row: segregation measured by the cumulative radial distribution of 
stars in different mass bins, measured from the position of the star in the region of highest stellar density. The lines, from light grey to black, show stars in the 
mass bins Af/M Q < 0.1, 0.1 < M/M s < 0.3, 0.3 < A//M Q < 1.0, 1.0 < M/Mq. Bottom row: segregation measured by the effective surface density 
of the most massive n stars (see the text for further explanation). The gray line is the median value for the surface density of n random stars, and the shaded 
regions show the ±lcr, ±2ct, and ±3cr values. The blue line is the surface density of the n most massive stars. 



Quantity / Distance range 


< 1 pc 


1 - 3pc 


3 - 10 pc 


10 - 30 pc 


30 - 100 pc 


> 100 pc 


Median mass [Mq] 


0.046 


0.033 


0.035 


0.047 


0.046 


0.050 


Upper quartile mass [Mq] 


0.14 


0.09 


0.11 


0.14 


0.15 


0.21 


Maximum mass [Mq] 


2.9 


5.3 


2.3 


1.7 


3.7 


3.7 


Velocity dispersion [km/s] 


0.8 


1.3 


0.7 


2.1 


5.7 


17 


Number objects 


276 


236 


223 


276 


169 


73 


Number binaries 


39 


20 


27 


32 


10 


4 


Binary fraction 


0.16 


0.09 


0.14 


0.13 


0.06 


0.06 



Table 2. Radial properties of the remnant cluster (radius 2 pc) and stellar halo (extending beyond 200 pc) at the end of the unperturbed W-body evolution 
with instantaneous gas removal (age 10 Myr). This table can be compared with Table 2 of Bate (2009a) for the properties of the cluster at the end of the 
hydrodynamical evolution. There is no evidence for radial mass segregation of the cluster in terms of the median, upper quartile, or maximum masses. Indeed, 
we note that some of the most massive stars in the cluster are present in each of the radial bins. The binary fraction is more or less uniform, except beyond 
30 pc where it drops by about a factor of two. The velocities of the ejected stars in the halo increase roughly linearly with distance beyond 3 pc (i.e. faster 
moving stars have travelled further). 



the timescale on which mass segregation manifests itself. This is 
the classical form of mass segregation, due to two-body effects in 
an almost spherical cluster, in contrast to the primordial segregation 
that is an imprint of star formation processes and clustered star for- 
mation. As such it should appear in cumulative distribution plots as 
a general trend, with more massive populations progressively more 
centrally concentrated. 

In Figure|9] we show the cumulative radial distributions for the 
main calculation and the three cases with gas potentials at 10 Myr. 
The case with Tb.i = 0.3 Myr shows the expected signal, with 
the distribution of the most massive stars beginning to be separated 
from the lower mass stars. The calculations with Tq.i = 1 Myr and 
without gas removal show segregation of the second-most massive 
bin as well. We stress that while the earliest (earlier than 1 Myr for 



these clusters) observations of segregation are a result of a form- 
ing cluster's morphology, this later type of mass segregation is ex- 
pected due to the dynamical maturity of these clusters, and does not 
contain information about the star formation process. 

3.3 The evolution of multiple systems 

|Bate|p009a| l quantified the fraction of stars and brown dwarfs that 
were in multiple systems using the multiplicity fraction, defined as 

B + T + Q 
S + B + T + Q' 

where S is the number of single stars within a given mass range 
and, B, T, and Q are the numbers of binary, triple, and quadru- 
ple systems, respectively, for which the primary has a mass in the 
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Figure 9. The dependence of the mass segregation at 10 Myr on the gas removal timescale. From left to right, the four calculations are: instantaneous gas 
removal, To.i = 0.3 Myr, Tq.i = 1 Myr, and no gas removal. The segregation is measured by the cumulative radial distribution of stars in different mass bins, 
measured from the position of the star in the region of highest stellar density. The lines, from light grey to black, show stars in the mass bins M/Mq< 0.1, 
0.1 ^ M/Mq < 0.3, 0.3 ^ M/M0< 1.0, 1.0 ^ M /Mq. The slow and no gas removal cases display classical mass segregation with more massive stars 
being progressively more centrally condensed (particularly in the two highest mass bins). 



same mass range. This differs from the companion star fraction, 
csf, that is also often used and where the numerator has the form 
.B+2r+3Q.|Bate|ch ose the multiplicity fraction following Hubber 
& Whitworth ( 2005 1 who pointed out that this measure is more ro- 
bust observationally in the sense that if a new member of a multiple 
system is found (e.g. a binary is found to be a triple) the quantity 
remains unchanged. It is also more robust for simulations too in 
the sense that if a high-order system decays because it is unstable 
the numerator only changes if a quadruple decays into two bina- 
ries (which is quite rare). Furthermore, if the denominator is much 
larger than the numerator (e.g. for brown dwarfs where the multi- 
plicity fraction is low) the production of a few single objects does 
not result in a large change to the value of MF. This was useful 
for Bate (2009al because many of the high-order systems in exis- 
tence at the end of the calculations were likely to undergo further 
dynamical evolution. Indeed, the long-term evolution of the mul- 
tiple systems produced by the main calculation of Bate (2009a) is 
the topic of this section. 

In Figure^] we plot the multiplicity fraction of the stars and 
brown dwarfs as a function of stellar mass for the calculation at the 
end of the hydrodynamical calculation (i.e. the beginning of the N- 
body evolution) and at the end of the main A'-body calculation at 
an age of 10 Myr (top left and top right panels, respectively). The 
top left panel is a reproduction of Figure 15 of|Bate (2009a}. The 
mass bins were chosen for comparison with the results of various 
observational surveys. The filled blue squares give the multiplicity 
fractions from the calculation while the surrounding blue hatched 
regions give the range in stellar masses over which the fraction is 
calculated and the la (68%) uncertainty on the multiplicity frac- 
tion. The black open boxes and their associated error bars and/or 
upper/lower limits give the results from a variety of observational 
surveys (see the figure caption). 

There is very little evolution of the multiplicity fraction dur- 
ing the main A'-body evolution of the cluster with instantaneous 
gas removal. Both the results at the end of the hydrodynamical cal- 
culation and those at the end of the A'-body evolution are in good 
quantitative and qualitative agreement with the observed multiplic- 
ity fraction (see Bate (2009al for a full discussion), which is found 
to be a strongly increasing function of primary mass. We also note 
that the decrease in binary fraction in the outer parts of the ejected 
halo of stars found by Bate (2009a) is maintained to ages of 10 Myr 
(Table[2](, but that the transition radius at which the binary fraction 
begins to decrease with distance has moved from 0. 15 pc to 30 pc 
as the halo has expanded. 

For the randomly-perturbed sets of initial conditions, although 



each calculation differs in detail with respect to which binaries 
and multiple systems survive or evolve dynamically, any varia- 
tions in the values of the multiplicity as a function of primary mass 
are smaller than the statistical uncertainties already plotted in Fig- 
ure [10] Typically, for each mass range one or two more or fewer 
multiple systems survives in each random realisation for primary 
masses < 0.50 Mq. For the more massive systems there tends to 
be somewhat more variation. For example, for the most massive 
systems (primary masses > 1.2 Mq) the multiplicity ranges from 
6/22 = 27% to 10/20 = 50%. But the statistical uncertainties 
are also largest for these multiplicities due to the small numbers of 
massive stars. 

The lower panels of Figure [10] give the multiplicity fractions 
for the two of the calculations aimed at investigating the depen- 
dence of the results on the gas removal time. The results from 
the To.r = 1 Myr (lower left panel) and no gas removal (lower 
right panel) calculations are given. It can be seen that the multiplic- 
ity fractions decrease slightly overall with increased gas removal 
timescales, but only marginally despite the fact that the remnant 
clusters become significantly more condensed and populous. Even 
in the extreme case of no gas removal, the decreases in the multi- 
plicity fractions for different primary masses between the instan- 
taneous gas removal and no gas removal cases are smaller than 
the statistical uncertainties. As for the randomly-perturbed sets of 
initial conditions, only the most massive systems (primary masses 
> 1.2 Mq) have a large change in the multiplicity fraction, but, 
again, the statistical uncertainties are also largest in this mass range 
due to the small numbers of massive stars. 



3.3.1 Triple and quadruple systems 

Although the multiplicity fraction does not evolve significantly, this 
is partially due to the fact that this measurement of the fraction 
of binaries and multiples was deliberately chosen to be relatively 
insensitive to the break up of dynamically-unstable multiple sys- 
tems. In fact, during the A'-body evolution of the unperturbed case 
with instantaneous gas removal, the number of quadruple systems 
drops from 25 to 8, while the number of triple systems remains the 
same (though the specific make up of the 23 systems changes). The 
number of binary systems increases slightly from 90 to 98 and the 
number of single objects also increases from 904 to 956 (Table[TJl. 
Much of this evolution occurs quite quickly, with the number of 
quadruple systems having dropped to 10 after 0.3 Myr of A'-body 
evolution (i.e. when the cluster age is 0.6 Myr). Thus, high-order 
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Figure 10. Multiplicity fraction as a function of primary mass. The top left panel gives the results at the end of the hydrodynamical calculation (0.3 Myr), 
while the other panels give the results at the end of the A'-body evolutions (10 Myr) for three cases: instantaneous gas removal (top right), To.i = 0.3 Myr 
(lower left), and with no gas removal (lower right). The blue filled squares surrounded by shaded regions give the results from the calculations with their 
statistical uncertainties. The open black squares with error bars and/or upper/lower limits give the observed multiplicity fractions from the surveys of Close et 
al. (2003), Basri & Reiners (2006), Fisher & Marcy (1992), Duquennoy & Mayor (1991), Preibisch et al. (1999) and Mason et al. (1998), from left to right. 
There is very little evolution of the multiplicity fraction over the 10 Myr of A'-body evolution, even for the extreme case of no gas removal. All sets of results 
are in reasonable agreement with the observed multiplicity as a function of primary mass. 



systems (particularly quadruple systems) are converted into lower- 
order systems through dynamical disruption, but the multiplicity 
fractions themselves are not significantly changed. 

The final overall frequencies of triple and quadruple systems 
are 2.1 and 0.7%, respectively. However, as discussed by |Bate| 
(2009a), these frequencies are a strong function of primary mass. 
At the end of the A-body evolution there are no quadruple systems 
with primary masses less than 0.2 Mq. There are two surviving 
very-low-mass (VLM) (0.03 — 0.10 Mq) triple systems out of a 
total of 331 systems (giving a frequency of 0.6%). At the end of 
the hydrodynamical calculation there were 3 triples and 2 quadru- 
ple systems with primaries in this mass range. Thus, the fraction 
of high-order VLM multiples has significantly decreased. For low- 
mass stars, there are 4 triple systems out of the 133 systems with 
primaries in the range 0.1 — 0.2 Mq (i.e. 3%) at the end of the A'- 
body evolution, whereas at the end of the hydrodynamical evolution 
there were 4 triples and 2 quadruples. The frequencies of triples or 
quadruples for M-dwarf primaries and masses 0.2 — 0.5 Mq is 9 
out of 96 systems (5 triples and 4 quadruples), very similar to the 7 
triples and 1 quadruple found from the 99 M-dwarf primaries sur- 
veyed by Fischer & Marcy ( 1992 1. At the end of the hydrodynami- 
cal evolution there were 5 triples and 10 quadruples with primaries 
in this mass range. For higher-mass primaries, the frequencies of 



triple and quadruple systems are 18% for 0.5 — 0.8 Mq and 29% 
for primary masses > 0.8 Mq, but with large statistical uncertain- 
ties. The number of high-order systems with these higher mass pri- 
maries is similar after the A-body evolution to that at the end of 
the hydrodynamical evolution, but whereas the number of triples 
has increased from 10 to 12, the number of quadruples has dropped 
from 1 1 to 4. 



The destruction of high-order multiple systems is qualitatively 
the same for the randomly-pertubed realisations and even for the 
calculations that include a gas potential. In Table[T] we give the total 
numbers of single, binary, triple and quadruple systems at the end 
of the hydrodynamical calculation and at the end of the four A-body 
calculations with different gas removal timescales. For gas removal 
timescales of 1 Myr or less, the numbers are statistically indistin- 
guishable. For the calculation without gas removal, the number of 
triples is lower, but on the other hand the number of quadruples is 
somewhat higher than the in the cases with gas removal. Closer in- 
spection (see the next section) shows that many of these quadruples 
are very wide systems in the halo. 
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Figure 11. The distributions of separations (semi-major axes) of multiple 
systems with stellar (left) and VLM (right) primaries. From top to bottom, 
the four rows give the results: at the end of the hydrodynamical calcula- 
tion (top row; 0.3 Myr; reproduced from Bate 2009a) and at the end of the 
/V-body evolution (10 Myr) with instantaneous gas removal (second row), 
gas removal timescale Trj.i = 1 Myr (third row), and no gas removal (bot- 
tom row). The To i = 0.3 Myr case has been omitted as the results lie 
between those of the instantaneous gas removal and Trj.i = 1 Myr cases. 
The solid, double hashed, and single hashed histograms give the orbital sep- 
arations of binaries, triples, and quadruples, respectively (each triple con- 
tributes two separations, each quadruple contributes three separations). In 
the stellar graph, the curve gives the G-dwarf separation distribution (scaled 
to match the area under the histogram in each case) from Duquennoy & 
Mayor (1991). In the VLM systems graph, the open black histogram gives 
the separation distribution (scaled to match the number in the 10-100 AU 
range of the calculated histogram in each case) of the known VLM multiple 
systems maintained by N. Siegler at http://vlmbinaries.org/ (last updated on 
February 4th, 2008). The vertical dotted line gives the resolution limit of 
the hydrodynamical calculation as determined by the gravitational soften- 
ing and accretion radii of the sink particles. 



3.3.2 Separation distributions 

Along with the evolution of the frequencies of multiple systems, we 
wish to understand how the characteristics of the multiple systems 
evolve. In Figure [TT] we present the separation (semi-major axis) 
distributions of the stellar (primary masses greater than 0.10 M©; 
left panels) and VLM multiples (right panels). These distributions 
are compared with the survey of Duquennoy & Mayor ( 199TJ and 
the listing of VLM multiples maintained by N. Siegler, C. Gelino, 
and A. Burgasser at http://vlmbinaries.org/ (last updated July 28, 
2009), respectively. The filled histograms give the separations of bi- 
nary systems, while the double hashed region adds the separations 
from triple systems (two separations for each triple, determined by 
sub-dividing the triple into a binary with a wider companion), and 
the single hashed region includes the separations of quadruple sys- 
tems (three separations for each quadruple which may be composed 
of two binary components or a triple with a wider companion). The 
top row gives the distributions at the end of the hydrodynamical 
calculation, while the second row gives the distributions at the end 
of the /V-body calculation (10 Myr) with instantaneous gas removal. 

We note that the gravitational softening and finite sink parti- 
cle accretion radii used in the hydrodynamical calculation of Bate 
(2009a) result in a pile up of stellar binaries with separations of a 
few AU. This was investigated by Bate] who re-ran part of the hy- 
drodynamical calculation without gravitational softening and using 
smaller accretion radii and, as expected, the pile up disappeared 
and a bell-shaped separation distribution was recovered (but with 
a much smaller number of objects because the re-run calculation 
could not be followed for as long). 

Comparing the upper two left-hand panels in Figure[TT|for the 
stellar multiple systems, we find that, as mentioned earlier, many 
of the quadruple systems are broken up, resulting in a marked de- 
crease in the contribution of the quadruple systems to the separa- 
tion distributions. There is also an increase in the number of triple 
systems with component separations < 10 AU. The other signif- 
icant change is in the number of very wide systems (separations 
10 4 — 10 5 AU). While at the end of the hydrodynamical calcula- 
tion there were only two separations in this range (one binary, and 
one from the wide component of a quadruple), at 10 Myr there are 7 
binaries, 5 triples, and 2 quadruples with separations in this range. 
These wide systems are found in the halo of ejected objects and 
are formed when two ejected objects have similar velocities and as 
they depart from the cluster, they happen to be weakly mutually 
bound. The overall effect of the break up of stellar quadruples with 
intermediate separations (1 — 1000 AU) and the formation of very 
wide systems is to significantly broaden the separation distribution. 
In fact, the separation distribution from 10 — 10 J AU appears to 
be even flatter than that observed by Duquennoy & Mayor ( 1991| >, 
and more consistent with the flat separation distribution of wide 
(500 — 5000 AU) binaries found for pre-main-sequence stars by 
|Kraus & Hillenbrand|p009) . 

For the VLM systems (primary masses < 0.1 Mq), the evo- 
lution is less dramatic (upper two right-hand panels of Figure [TTjl. 
Many of the triples and all of the quadruple systems present at the 
end of the hydrodynamical calculation have been broken up by the 
end of the /V-body evolution. However, these comprise a smaller 
fraction of the initial separation distribution than for the stellar sys- 
tems. There are two main effects of the long-term evolution. First, 
all the systems with separations greater than 60 AU at the end of 
the hydrodynamical calculation have been destroyed. Second, two 
extremely wide VLM binaries have been formed in the dispersing 
halo population: a 26,000 AU system consisting of 40 and 14 Mj 
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brown dwarfs and a 34,000 AU system consisting of 39 and 10 Mj 
brown dwarfs. This is good agreement with the observed separation 
distribution of VLM binaries: field VLM binaries are almost always 
found to have projected separatio ns ^ 40 AU {Martin et al.|2000| 
|Reid et al.|200"Tj|Close et al.|2003||Burgasser et al.|2007| > with only 
a handful of wider field systems found to date (220 AU; 1800 AU, 
5100 AU, 6 700 AU;|Billeres et al.|2005||Caballero|2007||Artigau| 
|et al.|2007| |Radigan et al.||2009) , while 4 out of 5 VLM binaries 
with separations in the 100 — 1000 AU range have been found in 
star- forming regi ons (|Luhm an 2004 ; Ja yawardhana & Ivanov|2006) 
|Close et al.|2007||Bejar et al.|2008> ' 

The evolution of the separation distributions discussed above 
is qualitatively similar for all the randomly-perturbed realisations, 
though quantitatively there is some variation. The number of wide 
(10 4 — 10 AU) systems displays significant variation. In the un- 
perturbed case, 7 binary, 5 triple and 2 quadruple separations fell 
within this separation range. The mean values for all 1 1 realisations 
are: 6.8 ± 3.5, 3.6 ± 2.7, and 4.5 ± 2.5 with the total number of 
separations in this range being 15 ± 7.2 with the minimum being 5 
and the maximum being 32. Thus, the unperturbed case is similar 
to the mean case, but sometimes there can be a lot more or a lot 
fewer wide multiple systems formed. This is not unexpected, since 
these systems require stars or systems to be ejected with similar 
velocities, something that is very sensitive to small perturbations. 

In the lower two rows of Figure [TT] we present results from 
the ,/V-body calculations with gas potentials at 10 Myr. The third 
row presents the results from the calculation with a gas removal 
timescale of Tb.i = 1 Myr, while the bottom row presents the re- 
sults from the calculation with no gas removal. The To.i = 0.3 Myr 
case has been omitted because, as might be expected, the results 
are nicely bracketed by the cases of instantaneous gas removal and 
Tb.i = 1 Myr. As expected, longer gas removal timescales destroy 
more wide systems (100 — 10 4 AU), with triple and quadruple sys- 
tems preferentially disrupted. However, there is very little differ- 
ence between the calculations with different gas removal timescales 
at separations less than 100 AU and the formation of significant 
numbers of very wide systems (separations > 10 4 AU) discussed 
above still occurs even without gas removal. For the VLM systems, 
there is also very little dependence on the gas removal timescale. 
In all cases, the main evolution is that the number of systems with 
separations of 100 — 10 4 AU decreases substantially and, in all but 
the no gas removal case, a few very wide (> 10 4 AU) systems are 
formed in the halo. Even in the extreme case of no gas removal 
at 10 Myr, the separation distribution is in good agreement with 
observations. 

3.3.3 Mass ratio distributions 

As we have seen, subsequent dynamical evolution of the |Bate| 
(2009a) cluster results in the formation of very wide systems and 
evolution of the separation of VLM binaries that is in good agree- 
ment with that observed. How do the mass ratio distributions 
evolve? 

In Figure [12] we present the mass ratio distributions of bi- 
naries with primaries that have masses ^ 0.5 Mq (left panels), 
are M-dwarfs with masses 0.1 ^ M < 0.5 Mq (centre panels), 
or are VLM objects (right panels). Again, the top panels are at 
the end of the hydrodynamical calculation (see also Figure 19 in 
|Bate| J2009a| Q, while the other panels give the results at the end 
of the iV-body evolutions at an age of 10 Myr. The case with in- 
stantaneous gas removal is given by the middle panels, while the 
lower panels give the results with no gas removal. We compare 



the M-dwarf mass ratio distribution to that of Fisch er & Marcy| 
{1992}, and the higher mass stars to the mass ratio distribution 
of solar-type stars obtained by Duquennoy & Mayor ( 1991 1. The 
VLM mass ratio distribution is compared with the listing of VLM 
multiples maintained by N. Siegler, C. Gelino, and A. Burgasser 
at http://vlmbinaries.org/ (last updated July 28, 2009). The figure 
only includes the mass ratios of binaries, but we include binaries 
that are components of triple and quadruple systems. A triple sys- 
tem composed of a binary with a wider companion contributes the 
mass ratio from the binary, as does a quadruple composed of a triple 
with a wider companion. A quadruple composed of two binaries or- 
biting each other contributes two mass ratios - one from each of the 
binaries. 

As noted by Bate at the end of the hydrodynamical calcula- 
tion, the ratio of near-equal mass systems to systems with dissimilar 
masses decreases going from VLM objects to M dwarfs in a similar 
way to the observed mass ratio distributions, but that the trend is not 
as strong as in the observed systems. Specifically, 71% of the VLM 
binaries had M2/M1 > 0.6 while for primary masses 0.1 — 0.5 
Mq the fraction was only 51%. The M-dwarf mass ratio distribu- 
tion is consistent with Fischer & Marcyj s distribution. The VLM 
binaries, although biased towards equal-mass systems, are not as 
strongly biased as is observed. The main problem with the hydro- 
dynamical binary mass ratio distributions, however, is that amongst 
'solar-type' binaries with primary masses greater than 0.5 Mq the 
proportion with mass ratios M2/M1 < 0.5 is much lower than that 
found by Duquennoy & Mayor) 

At the end of the ,/V-body evolution with instantaneous gas re- 
moval, the mass ratio distributions of the VLM and M-dwarf bina- 
ries have not changed significantly. 75% of the VLM binaries have 
M2/M1 > 0.6 and 55% of the M-dwarfs. However, there is a sub- 
stantial change in the mass ratio distribution of solar-type binaries. 
The number of near-equal mass binaries has decreased somewhat, 
while the number of unequal mass binaries has substantially in- 
creased. Whereas only 9/34 = 26% had M 2 /M± < 0.4 at the end 
of the hydrodynamical calculation, the fraction is 15/36 = 42% at 
the end of the A'-body evolution. 

As with the evolution of the separation distributions, the evo- 
lution of the mass ratio distributions discussed above is qualita- 
tively similar for all the randomly-perturbed realisations. The num- 
ber of unequal-mass solar-type binaries displays some variation. In 
the unperturbed case there were 15 with M2/M1 < 0.4 (out of 
36), while the mean value is slightly lower at 10.7 ± 2.3, and the 
minimum and maximum values are 7 and 15, respectively. There is 
also some variation in the number of unequal-mass VLM binaries. 
The unperturbed case left 9 with M2/M1 < 0.6 (out of 36), while 
the mean value is 8.7 ± 2.0 and the minimum and maximum values 
are 6 and 12, respectively. 

The mass ratio distributions resulting from the TV-body calcu- 
lations with gas potentials also differ very little from those obtained 
with instantaneous gas removal. Even the extreme case of no gas 
removal differs only slightly from the instantaneous case (compare 
the middle and lower rows of Figure |T2). 

To illustrate what has caused the increase in unequal-mass 
solar-type binaries in the case with instantaneous gas removal, 
we discuss the formation of the three binaries with mass ratios 
0.2 < M2/M1 < 0.4. One of these systems was formed when 
a triple consisting of a 2.5 AU binary containing 0.96 and 0.73 Mq 
stars with an 0.21 Mq companion at 12 AU broke up to leaving the 
0.73 and the 0.21 Mq stars in a tight (0.1 AU separation) binary 
with a low-mass ratio. Another was formed when a quadruple sys- 
tem consisting of two binaries with 21 AU and 5 AU semi-majors 
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Figure 12. The mass ratio distributions of binary systems with stellar primaries in the mass ranges M\ > 0.5 Mq (left) and Mi = 0.1 — 0.5 Mq 
(centre) and VLM primaries (right; Mi < 0.1 Mq). The top panels give the results at the end of the hydrodynamical calculation (0.3 Myr; Bate 2009), 
while the results at the end of the iV-body evolutions (10 Myr) with instantaneous gas removal and no gas removal are given in the middle and lower panels, 
respectively. The solid black lines give the observed mass ratio distributions of Duquennoy & Mayor (1991) for G dwarfs (left), Fischer & Marcy (1992) for 
Mi = 0.3 — 0.57 Mq (centre, solid line) and Mi = 0.2 — 0.57 Mq (centre, dashed line), and of the known VLM binary systems maintained by N. Siegler, 
C. Gelino, and A. Burgasser at http://vlmbinaries.org/ (right). The observed mass ratio distributions have been scaled so that the areas under the distributions 
(M2/M1 = 0.4 — 1.0 only for the centre panels) match those from the simulations. In all cases, VLM binaries are biased towards equal masses when 
compared with M dwarf binaries (primary masses in the range Mi = 0.1 — 0.5 Mq). At 10 Myr with instantaneous gas removal, 75% of the VLM binaries 
have Mq, /M\ > 0.6 while for the M dwarf binaries the fraction is only 55%. There is only weak evolution of the VLM and M-dwarf mass ratio distributions 
during the /V-body evolution, even in the most extreme case of no gas removal. However, with fast gas removal the /V-body evolution results in a significant 
increase in the fraction of unequal mass solar-type binaries (the fraction of systems with mass ratios M2/M1 < 0.4 increases from 26% to 42% in the case 
with instantaneous gas removal). 



axes orbiting each other with a 95-AU semi-major axis broke up 
forming a new binary with a mass ratio of M2/M1 = 0.31 con- 
sisting of one component from each of the original binaries in a 12- 
AU semi-major axis. The last was formed when 1.74 and 0.40 Mq 
stars that had been ejected from the cluster with similar velocities 
happened to find themselves weakly mutually bound in a 4100-AU 
binary. As was seen in the previous section, a significant number 
of wide systems are formed in this manner and they tend to con- 
tain (although not exclusively) of one of the more massive stars 
(> 0.5 Mq), as might be expected for a wide pair of objects to 



be bound. Another low-mass ratio example consists of 0.70 and 
0.10 M Q stars in a 7500- AU binary. 

The formation of low-mass ratio solar-type binaries during the 
V-body evolution with instantaneous gas removal goes some way 
to addressing the problem that hydrodynamical simulations of star 
formation appear to under-produce unequal-mass solar-type bina- 
ries (|Delgado-Donate et al.|2004[[Bate|2009a[parke|2009[ l. [Batel 
( 2009a ) suggested that if the mass ratios of triples and quadruples 
were included this might also help to address the disagreement (i.e. 
essentially arguing that some of Duquennoy & Mayor s unequal 
mass wide binaries might have in fact been unidentified triples or 
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Figure 13. The mass ratio distributions of binary, triple and quadruple sys- 
tems with stellar primaries with masses Mi > 0.5 Mq at the end (10 Myr) 
of the unperturbed /V-body evolution with instantaneous gas removal, ex- 
cluding VLM companions (masses < 0.1 Mq. The double hashed region 
gives the mass ratio distribution for binaries only. The single hashed region 
includes the mass ratios of triple and quadruple systems taking the mass 
ratio of the two most widely separated components. The solid black line 
gives the fit to the observed mass ratio distribution of Duquennoy & Mayor 
(1991) for G dwarfs. The observed mass ratio distribution has been scaled 
so that the areas under the distribution matches that from the simulation. 



quadruples). In Figure [13] we plot the mass ratio distribution at 
the end of the A'-body evolution of systems with primary masses 
> 0.5 Mq including the mass ratios between the wide compo- 
nents of triples and quadruple systems, but excluding VLM com- 
panions. For example, for a triple, the mass ratio is that between 
the close binary and the wider companion. For a quadruple consist- 
ing of two binaries, we include the mass ratio of the two binaries. 
We exclude VLM companions because the survey of Duquennoy 
|& Mayor] ( |1991) was not sensitive to VLM objects. The inclusion 
of triple and quadruple systems, in addition to the formation of 
unequal-mass binaries during the A'-body evolution, does not com- 
pletely remove the discrepancy between the simulation and the sur- 
vey of Duquennoy & Mayor ( 1991 ), but it does go a long way to- 
wards addressing it. We also note that the recent surveys of binaries 
with separations 5 — 5000 AU in star-forming regions ( Kraus et al. 
|2008[|Kraus & H illenbrand 2009 1 find uniform mass ratio distribu- 
tions for primary masses ^ 0.75 Mq, unlike that of Duquennoy &| 



|Mayor| ( (l99l) . However, [Kraus & Hillenbrand| ( |2009| > also found 
an indication of a variation between different star-forming regions 
with 9/12 of the wide binaries (500 — 5000 AU) in Taurus hav- 
ing M2/M1 > 0.75 while in Upper Sco 6/11 of the binaries had 
M2/M1 < 0.25. In summary, the mass ratio distribution of bina- 
ries warrants further investigation, both observationally and theo- 
retically. 



4 DISCUSSION 

For the first time, we have studied the evolution following gas dis- 
persal of a stellar cluster produced from a hydrodynamical calcu- 
lation of the collapse and fragmentation of a turbulent molecular 
cloud. Since only one hydrodynamical star formation calculation 
that forms a substantial cluster (2> 100 stars) and resolves binary 
and multiple stellar systems has ever been performed (Bate 2009a) 
our results are, by necessity, limited to one particular type of clus- 
ter in terms of total mass, stellar density, etc. However, while it is 
difficult to extrapolate our results more widely (e.g. to more or less 
massive clusters), the calculations presented here at least illustrate 



the types of evolution that are likely to occur in other types of clus- 
ters. 



4.1 Cluster evolution during gas dispersal 

The cluster studied here results from the hydrodynamical collapse 
and fragmentation of a dense turbulent molecular cloud and, just 
before gas dispersal, it is very compact (half-mass radius 0.05 pc), 
with a very high central stellar density (~ 10 6 pc 3 ), surrounded by 
a halo of ejected objects. Despite this initial state, upon removal of 
the gas (which makes up 62% of the total mass), a bound remnant 
cluster, that contains about 30-40% of the mass and number of stars, 
expands to a radius of ~ 1 — 2 pc over a period of ~ 4 — 10 Myr 
(depending on the gas removal timescale). This evolution is consis- 
tent with the many past A'-body studies of the effects gas explusion 
on cluster evolution which have started from idealised initial condi- 
tions. For example, Baumgardt & Kroupa ( 2007 1 performed a wide 
variety of A'-body simulations of the cluster dissolution following 
gas dispersal with varing star formation efficiencies. They found 
that their surviving clusters expanded by up to a factor of 10 in size. 
The obvious question raised by this sort of evolution is whether or 
not real stellar clusters may originate from such compact configu- 
rations? 

The amount of gas removed from the cluster we study here is 
large and only 30 — 40% of the total number of stars remain in the 
remnant cluster. However, these numbers are not unreasonable. A 
star formation efficiency of ~ 40% has been inferred for the ONC 
^Hillenbrand & Hartmann|1998| >, while |Weidner et al.|p007) find 
observationally that low-mass clusters keep at most 50% of their 
stars during gas loss. 

Recently, the observational evidence for significant expansion 
of young clusters has also been mounting. |Scally et al.| p005[ > 
showed that if the current estimates of the virial ratio of the ONC 
(which show that the cluster is unbound) are correct, the size and 
age of the ONC imply either that it became unbound only very re- 
cently, or else that it has expanded quasi-statically. In the latter case, 
they infer that its initial central density may have exceeded its cur- 



rent value of ~ 10 4 pc" 3 ( |McCaughrean & Stauffer|l994 



Hillen- 



|brand & Hartmann 1998 ) by 1-2 orders of magnitude which would 
make it similar to the cluster studied here. Parker et al. ( 20091 also 
argue that the central region of the ONC was 2 orders of magni- 
tude denser in the past with a half-mass radius of only 0.1-0.2 pc 
based on its current population of wide binaries. Moving to clus- 
ters other than the ONC, Bastian et al. ( 2008 1 studied several dozen 
Galactic and extragalactic stellar clusters and found evidence for 
rapid expansion of the cluster cores during the first 20 Myr of their 
evolution. Examining the evolution of cluster sizes with age, they 
found that clusters may begin with core radii ^ 0.1 pc in size. 

How else might such young cluster evolution be tested? From 
their A'-body simulations, Baumgardt & Kroupa (20071 found 
that if gas explusion is rapid, the stars acquire strongly radially 
anisotropic velocity dispersions outside their half-mass radii. Simi- 
larly, in the cluster evolution presented here, the ejected halo popu- 
lation which comprises > 60% of the stars and stellar mass are on 
unbound essentially radial trajectories with velocities that increase 
with distance from the cluster remnant (faster moving stars have 
travelled further). 

Finally, as noted in Section [3.1.2| a common feature of the A'- 
body calculations is for the most massive binary to be ejected from 
the remnant cluster, sometimes accompanied by a small group of ~ 
20 companion stars. Many Herbig Ae/Be stars are associated with 
small groups of stars ( Hillenbra ndp994] [T995] Hillenbra nTeTaT] 
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|1995[ >. Since the ejection of such a group occurs in the several of 
our randomly-perturbed realizations, we suggest that this might be 
a significant formation mechanism of Herbig Ae/Be stars in small 
groups as opposed to an intermediate-mass star forming with a few 
other stars in an isolated molecular cloud. 



that is usually present at the end of the hydrodynamical star forma- 
tion phase jDelgado-Donate et al.|2004{|Bate|2009a||Clarke|2009| ). 
Thus, in summary, the formation of significant numbers of wide 
multiple systems during cluster dissolution potentially solves sev- 
eral problems simultaneously. 



4.2 The evolution of binary and multiple stellar systems 

As noted in Section[T| many previous A'-body simulations have in- 
vestigated the evolution of binary properties, particularly binary 
frequency and separation distributions, during the early evolution 
of a young cluster (e.g |Kroupa|1995a|b|[T99"8l|Parker et al.|2009| >. 
Many of the conclusions reached from these studies are also appli- 
cable to the cluster we have studied here. In particular, we find a 
cluster remnant surrounded by a halo of ejected stars with a lower 
binary fraction ( Kroupa 1995b). We also find that some of the most 
massive stars can be ejected and reach large distances from the clus- 
ter within 10 My r (|Kroupa|1 998 ). However, unlike [Kroupa| \ 1 995a) 
and |Parker et al. l (2009l > we do not find that many primordial binary 
systems are disrupted. Rather, the vast majority of binary systems 
formed during the star formation process (i.e. the hydrodynamical 
evolution) survive the gas dispersal phase with little evolution of bi- 
nary frequency or their separation distribution. This is because stars 
do not form in the hydrodynamical simulations with 100% of them 
being in binaries, nor do they form with a flat distribution in log- 
separation. Rather, star formation and dynamical interaction cannot 
be separated - they proceed simultaneously, resulting in a 'primor- 
dial' binary population that is already somewhat stable against fur- 
ther evolution. To try to separate star formation and the properties 
of primordial binary properties is likely to severely over-estimate 
the importance of the dynamical evolution of the binary population 
early in the life of a stellar cluster. 

By contrast, we note that higher-order multiple systems 
(which are usually not included in purely A'-body studies at all) do 
undergo significant evolution during the early evolution of a stellar 
cluster. In particular, many quadruple systems that result from the 
star formation process are found to decay into binaries or triples. 
Thus, we assert that observational studies of higher-order multi- 
ple systems in very young stellar clusters are likely to find a much 
higher proportion of high-order multiples than is found in open 
clusters or the field population (e.g. Reipurth 200T)]>. 

Finally, as presented in Section 3.3 we find that a significant 
fraction of very wide (separations 10 — 10 5 AU) multiple systems 
(around half of which are actually triples or quadruples) are pro- 
duced in the halo of objects ejected from the cluster. This has also 
been noted in recent pure A'-body simulations of dispersing clusters 
( |Kouwenhoven et al.|20 09 1. Thus, even though the hydrodynami- 
cal star formation process in a very dense cluster such as that dis- 
cussed here (half-mass radius only 0.05 pc, or 10 4 AU at the end 
of the star formation) cannot not produce very wide systems by 
direct fragmentation, such wide systems can be produced during 
cluster dissolution in numbers that are large enough to explain the 
observed field population of wide systems (Duquennoy & Mayor 
|1991| >. Such formation of wide systems during cluster dissolution 
potentially solves the problem i Parker et al.|2009| of how a broad 
log-separation distribution (e.g. Duquennoy & Mayor|19"9T l can be 
obtained if most stars form in clusters jAdams & M yers 200l| |Lada| 
|& Lada|2003| l which cannot contain wide binaries. Moreover, many 
of these wide binary systems and the binaries formed from the de- 
cay of higher-order multiple systems have unequal mass compo- 
nents. This also helps to fix the deficit of unequal-mass solar-type 
binaries (relative to the observations of Duquennoy & Mayor 1991 1 



5 CONCLUSIONS 

We have presented the results of A'-body calculations that take the 
end point of the hydrodynamical star cluster formation calculation 
of Bate (2009a I and evolve the stellar cluster to an age of 10 Myr. 
The calculations are unique in the sense that they begin from a star 
cluster that has self-consistently formed from the hydrodynamical 
collapse and fragmentation of a molecular cloud to form a stellar 
cluster containing binary and higher-order multiple systems rather 
than beginning from arbitrary initial conditions. Although we are 
limited to studying one particular type of cluster in terms of the 
number of stars, stellar density, etc, we find many of the same evo- 
lutionary processes that have been studied in past A'-body calcula- 
tions that begin from arbitrary initial conditions. For example, we 
find that a bound remnant cluster containing « 30 — 40% of the 
stars by number or mass remains despite the low star formation ef- 
ficiency (38%) and this cluster is surrounded by a halo of ejected 
stars. The remnant cluster expands from an initial radius of less than 
0.05 pc to « 1 — 2 pc over a period of w 4 — 10 Myr (depending on 
the gas removal timescale). This result, along with other recent nu- 
merical and observational studies, suggest that young clusters may 
begin from much higher stellar densities (10 5 — 10 7 pc~ J ) than 
usually assumed. 

When investigating mass segregation in our young clusters, 
we differentiate between primordial mass segregation, which is a 
result of the star formation process, and classical mass segrega- 
tion, which is due to dynamical relaxation. We find that primordial 
mass segregation is only evident for the few most massive stars (in 
agreement with observations of the Orion Nebula Cluster), but only 
lasts a short time (~ 1 Myr). Later, classical mass segregation may 
develop with higher-mass stars being increasingly centrally con- 
densed, but this depends on the gas removal timescale. Slower gas 
removal results in shorter relaxation times and, thus, more classical 
mass segregation. 

Although many of the global cluster evolutionary phases have 
been discussed before in the literature, we also find some new evo- 
lutionary processes at work. We find that in the majority of our cal- 
culations, some of the most massive stars in the cluster are ejected 
from the bound cluster and are accompanied by a small group of 
companion stars. Thus, we propose that this mechanism might be 
responsible for the formation of Herbig Ae/Be stars that are accom- 
panied by small stellar groups. 

Turning to the evolution of binary and multiple systems, we 
note that past A'-body studies of the evolution of young clusters 
that begin with high binary frequencies and wide separation distri- 
butions have emphasised the rapid destruction of binary systems. 
However, in our calculations the binary frequency remains almost 
unchanged throughout the A'-body evolution. We emphasise that 
the processes of star formation and the formation of multiple stel- 
lar systems cannot be separated during the formation of a clus- 
ter. Instead, they occur simultaneously, resulting in a 'primordial' 
binary population that is already somewhat stable against further 
evolution. That said, we do find that the higher-order multiple sys- 
tems produced by the star formation process, particularly quadru- 
ples, evolve significantly with most decaying into lower-order sys- 
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terns. This is potentially observable by comparing the frequencies 
of triple and quadruple systems in young star-forming regions with 
those found in open clusters. We also find very wide (10 4 -10 5 AU) 
binary and higher-order multiple systems are formed in the dispers- 
ing halo of ejected objects when objects happen to be expelled with 
similar velocities and find themselves weakly bound. These wide 
binaries are formed in sufficient numbers to explain the observed 
separation distribution of solar-type stars despite the fact that the 
original star-forming cloud is too small to form such systems by 
direct fragmentation. Finally, whereas the cluster at the end of the 
hydrodynamical evolution was deficient in unequal-mass solar-type 
binaries compared with observations, we find that many of the wide 
binaries formed in the halo and those binaries resulting from the 
decay of higher-order multiple systems have unequal mass compo- 
nents. Thus, these processes may help to bring the mass ratio dis- 
tributions of binaries produced by hydrodynamical simulations into 
agreement with the observed mass ratio distributions of solar-type 
stars. 
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